#ifndef GUARD_maths
#define GUARD_maths

#include <math.h>
const float eps=1e-6f;
const double PI=3.1415926535897932384626433832795;

float Eps(float l)
{
	float k=1.f;
	float m=float(fabs(l));
	if (m>1e6f) k=m/1e6f;
	else if (m<1e-3f) k=m*1e3f+1e-20f;
	return k*eps;
}
float Eps(double l)
{
	float k=1.f;
	float m=float(fabs(l));
	if (m>1e6f) k=m/1e6f;
	else if (m<1e-3f) k=m*1e3f+1e-20f;
	return k*eps;
}

struct Vector3;

struct point3
{
	float x,y,z;

	point3(float xx=0.,float yy=0.,float zz=0.):x(xx),y(yy),z(zz) {}
	point3(const point3& p):x(p.x),y(p.y),z(p.z) {}
	point3& operator= (const point3& p)
	{
		x=p.x;
		y=p.y;
		z=p.z;
		return *this;
	}
	
	point3 operator+ (const point3& p) const
	{
		return point3(x+p.x,y+p.y,z+p.z);
	}
	point3 operator- (const point3& p) const
	{
		return point3(x-p.x,y-p.y,z-p.z);
	}
	point3 operator+ (const Vector3& v) const;
	point3 operator- (const Vector3& v) const;
	point3 operator* (const float& n) const
	{
		return point3(x*n,y*n,z*n);
	}
	point3 operator/ (const float& n) const
	{
		return point3(x/n,y/n,z/n);
	}
	bool operator== (const point3& p) const
	{
		float E=Eps((fabs(x)+fabs(y)+fabs(z))/3);
		if (fabs(x-p.x)<=E && fabs(y-p.y)<=E && fabs(z-p.z)<=E) return true;
		else return false;
	}
/*	bool operator== (const point3& p) const
	{
		if (fabs(x-p.x)<=(x?fabs(x):fabs(p.x))*eps && fabs(y-p.y)<=(y?fabs(y):fabs(p.y))*eps &&
				fabs(z-p.z)<=(z?fabs(z):fabs(p.z))*eps) return true;
		else return false;
	}
	bool operator== (const point3& p) const
	{
		if (fabs(x-p.x)<=eps && fabs(y-p.y)<=eps && fabs(z-p.z)<=eps) return true;
		else return false;
	}
	bool operator!= (const point3& p) const
	{
		if (fabs(x-p.x)<eps && fabs(y-p.y)<eps && fabs(z-p.z)<eps) return false;
		else return true;
	}
*/
	bool operator!= (const point3& p) const
	{
		float E=Eps((fabs(x)+fabs(y)+fabs(z))/3);
		if (fabs(x-p.x)<=E && fabs(y-p.y)<=E && fabs(z-p.z)<=E) return false;
		else return true;
	}
};

struct point4
{
	float x,y,z,f;

	point4(float xx=0.,float yy=0.,float zz=0.,float ff=0.):x(xx),y(yy),z(zz),f(ff) {}
	point4(const point4& p):x(p.x),y(p.y),z(p.z),f(p.f) {}
	point4& operator= (const point4& p)
	{
		x=p.x;
		y=p.y;
		z=p.z;
		f=p.f;
		return *this;
	}
	
	point4 operator+ (const point4& p) const
	{
		return point4(x+p.x,y+p.y,z+p.z,f+p.f);
	}
	point4 operator- (const point4& p) const
	{
		return point4(x-p.x,y-p.y,z-p.z,f-p.f);
	}
	point4 operator* (const float& n) const
	{
		return point4(x*n,y*n,z*n,f*n);
	}
	point4 operator/ (const float& n) const
	{
		return point4(x/n,y/n,z/n,f/n);
	}
	bool operator== (const point4& p) const
	{
		float E=Eps((fabs(x)+fabs(y)+fabs(z)+fabs(f))/4);
		if (fabs(x-p.x)<=E && fabs(y-p.y)<=E && fabs(z-p.z)<=E && fabs(f-p.f)<=E) return true;
		else return false;
	}
/*	bool operator== (const point4& p) const
	{
		if (fabs(x-p.x)<=(x?fabs(x):fabs(p.x))*eps && fabs(y-p.y)<=(y?fabs(y):fabs(p.y))*eps &&
			fabs(z-p.z)<=(z?fabs(z):fabs(p.z))*eps && fabs(f-p.f)<=(f?fabs(f):fabs(p.f))*eps) return true;
		else return false;
	}
	bool operator== (const point4& p) const
	{
		if (fabs(x-p.x)<=eps && fabs(y-p.y)<=eps && fabs(z-p.z)<=eps && fabs(f-p.f)<=eps) return true;
		else return false;
	}
	bool operator!= (const point4& p) const
	{
		if (fabs(x-p.x)<=eps && fabs(y-p.y)<=eps && fabs(z-p.z)<=eps && fabs(f-p.f)<=eps) return false;
		else return true;
	}
*/
	bool operator!= (const point4& p) const
	{
		float E=Eps((fabs(x)+fabs(y)+fabs(z)+fabs(f))/4);
		if (fabs(x-p.x)<=E && fabs(y-p.y)<=E && fabs(z-p.z)<=E && fabs(f-p.f)<=E) return false;
		else return true;
	}
};

struct point2
{
	float x,y;

	point2(float xx=0.,float yy=0.):x(xx),y(yy) {}
	point2(const point2& p):x(p.x),y(p.y) {}
	point2& operator= (const point2& p)
	{
		x=p.x;
		y=p.y;
		return *this;
	}
	
	point2 operator+ (const point2& p) const
	{
		return point2(x+p.x,y+p.y);
	}
	point2 operator- (const point2& p) const
	{
		return point2(x-p.x,y-p.y);
	}
	point2 operator* (const float& n) const
	{
		return point2(x*n,y*n);
	}
	point2 operator/ (const float& n) const
	{
		return point2(x/n,y/n);
	}
	bool operator== (const point2& p) const
	{
		float E=Eps((fabs(x)+fabs(y))/2);
		if (fabs(x-p.x)<=E && fabs(y-p.y)<=E) return true;
		else return false;
	}
/*	bool operator== (const point2& p) const
	{
		if (fabs(x-p.x)<=(x?fabs(x):fabs(p.x))*eps && fabs(y-p.y)<=(y?fabs(y):fabs(p.y))*eps) return true;
		else return false;
	}
	bool operator== (const point2& p) const
	{
		if (fabs(x-p.x)<=eps && fabs(y-p.y)<=eps) return true;
		else return false;
	}
	bool operator!= (const point2& p) const
	{
		if (fabs(x-p.x)<eps && fabs(y-p.y)<eps) return false;
		else return true;
	}
*/
	bool operator!= (const point2& p) const
	{
		float E=Eps((fabs(x)+fabs(y))/2);
		if (fabs(x-p.x)<=E && fabs(y-p.y)<=E) return false;
		else return true;
	}
};

struct Vector2
{
	float x,y;

	Vector2(float xx=0.,float yy=0.):x(xx),y(yy) {}
	Vector2(const Vector2& v):x(v.x),y(v.y) {}
	Vector2(const point2& p):x(p.x),y(p.y) {}
	Vector2(const point2& p1,const point2& p2):x(p2.x-p1.x),y(p2.y-p1.y) {}
	Vector2& operator= (const Vector2& v)
	{
		x=v.x;
		y=v.y;
		return *this;
	}
	Vector2 operator+ (const Vector2& v) const
	{
		return Vector2(x+v.x,y+v.y);
	}
	Vector2& operator+= (const Vector2& v)
	{
		x+=v.x;
		y+=v.y;
		return *this;
	}
	Vector2 operator- (const Vector2& v) const
	{
		return Vector2(x-v.x,y-v.y);
	}
	Vector2& operator-= (const Vector2& v)
	{
		x-=v.x;
		y-=v.y;
		return *this;
	}
	Vector2 operator* (const float& n) const
	{
		return Vector2(x*n,y*n);
	}
	Vector2 operator/ (const float& n) const
	{
		return Vector2(x/n,y/n);
	}
	float operator* (const Vector2& v) const		//Scalar product
	{
		return x*v.x+y*v.y;
	}
};

//Normal to the line
Vector2 normal2d(const point2& p1, const point2& p2)
{
	return Vector2(p2.y-p1.y,p1.x-p2.x);
}
Vector2 normal2d(const Vector2& v)
{
	return Vector2(v.y,-v.x);
}

//Square of the distance between 2 points in 2D
float dist2sq(float x1,float y1,float x2,float y2)
{
	return ( (x2-x1) * (x2-x1) + (y2-y1) * (y2-y1) );
}

//Square of the distance between 2 points in 3D
float dist3sq(float x1,float y1,float z1,float x2,float y2,float z2)
{
	return ( (x2-x1) * (x2-x1) + (y2-y1) * (y2-y1) + (z2-z1) * (z2-z1) );
}
float dist3sq(const point3& p1,const point3& p2)
{
	return ( (p2.x-p1.x) * (p2.x-p1.x) + (p2.y-p1.y) * (p2.y-p1.y) + (p2.z-p1.z) * (p2.z-p1.z) );
}
float dist4sq(const point4& p1,const point4& p2)
{
	return ( (p2.x-p1.x) * (p2.x-p1.x) + (p2.y-p1.y) * (p2.y-p1.y) + (p2.z-p1.z) * (p2.z-p1.z) + (p2.f-p1.f) * (p2.f-p1.f) );
}

//The verification of (x,y) to belonging to a triangle
/*int p_in_tr(float x,float y,float x1,float y1,float x2,float y2,float x3,float y3)
{
	double a,b,c,l12,l13,l23,l01,l02,l03;

	a=( (y3-y1) * (x2-x1) - (y2-y1) * (x3-x1) ) * ( (y-y1) * (x2-x1) - (y2-y1) * (x-x1) );
	b=( (y1-y3) * (x2-x3) - (y2-y3) * (x1-x3) ) * ( (y-y3) * (x2-x3) - (y2-y3) * (x-x3) );
	c=( (y2-y3) * (x1-x3) - (y1-y3) * (x2-x3) ) * ( (y-y3) * (x1-x3) - (y1-y3) * (x-x3) );
	l01=sqrt(dist2sq(x1,y1,x,y));
	l02=sqrt(dist2sq(x2,y2,x,y));
	l03=sqrt(dist2sq(x3,y3,x,y));
	l12=sqrt(dist2sq(x1,y1,x2,y2));
	l23=sqrt(dist2sq(x2,y2,x3,y3));
	l13=sqrt(dist2sq(x1,y1,x3,y3));
	if ( (a>0.) && (b>0.) && (c>0.) || (l01+l02-l12<eps/1000) || (l01+l03-l13<eps/1000) || (l02+l03-l23<eps/1000) ) return 1;
	else return 0;
}*/
int p_in_tr(float x,float y,float x1,float y1,float x2,float y2,float x3,float y3)
{
	Vector2 v12=Vector2(x2-x1,y2-y1);
	Vector2 v13=Vector2(x3-x1,y3-y1);
	Vector2 v23=Vector2(x3-x2,y3-y2);
	Vector2 n12=normal2d(v12);
	Vector2 n13=normal2d(v13);
	Vector2 n23=normal2d(v23);
	Vector2 v1=Vector2(x-x1,y-y1);
	Vector2 v2=Vector2(x-x2,y-y2);

	float E=Eps((fabs(x)+fabs(y)+fabs(x1)+fabs(y1)+fabs(x2)+fabs(y2)+fabs(x3)+fabs(y3))/8);

	if ((n12*v13)*(n12*v1)<-E) return 0;
	if ((n13*v12)*(n13*v1)<-E) return 0;
	if ((n23*v12)*(n23*v2)>E) return 0;

	return 1;
}

struct Vector3
{
	float x,y,z;

	Vector3(float xx=0.,float yy=0.,float zz=0.):x(xx),y(yy),z(zz) {}
	Vector3(const Vector3& v):x(v.x),y(v.y),z(v.z) {}
	Vector3(const point3& p):x(p.x),y(p.y),z(p.z) {}
	Vector3(const point3& p1,const point3& p2):x(p2.x-p1.x),y(p2.y-p1.y),z(p2.z-p1.z) {}
	Vector3& operator= (const Vector3& v)
	{
		x=v.x;
		y=v.y;
		z=v.z;
		return *this;
	}
	Vector3 operator+ (const Vector3& v) const
	{
		return Vector3(x+v.x,y+v.y,z+v.z);
	}
	Vector3& operator+= (const Vector3& v)
	{
		x+=v.x;
		y+=v.y;
		z+=v.z;
		return *this;
	}
	Vector3 operator- (const Vector3& v) const
	{
		return Vector3(x-v.x,y-v.y,z-v.z);
	}
	Vector3& operator-= (const Vector3& v)
	{
		x-=v.x;
		y-=v.y;
		z-=v.z;
		return *this;
	}
	Vector3 operator* (const float& n) const
	{
		return Vector3(x*n,y*n,z*n);
	}
	Vector3 operator/ (const float& n) const
	{
		return Vector3(x/n,y/n,z/n);
	}
	float operator* (const Vector3& v) const		//Scalar product
	{
		return x*v.x+y*v.y+z*v.z;
	}
};

point3 point3::operator+ (const Vector3& v) const
{
	return point3(x+v.x,y+v.y,z+v.z);
}
point3 point3::operator- (const Vector3& v) const
{
	return point3(x-v.x,y-v.y,z-v.z);
}

struct Vector4
{
	float x,y,z,f;

	Vector4(float xx=0.,float yy=0.,float zz=0.,float ff=0):x(xx),y(yy),z(zz),f(ff) {}
	Vector4(const Vector4& v):x(v.x),y(v.y),z(v.z),f(v.f) {}
	Vector4(const point4& p):x(p.x),y(p.y),z(p.z),f(p.f) {}
	Vector4(const point4& p1,const point4& p2):x(p2.x-p1.x),y(p2.y-p1.y),z(p2.z-p1.z),f(p2.f-p1.f) {}
	Vector4& operator= (const Vector4& v)
	{
		x=v.x;
		y=v.y;
		z=v.z;
		f=v.f;
		return *this;
	}
	Vector4 operator+ (const Vector4& v) const
	{
		return Vector4(x+v.x,y+v.y,z+v.z,f+v.f);
	}
	Vector4& operator+= (const Vector4& v)
	{
		x+=v.x;
		y+=v.y;
		z+=v.z;
		f+=v.f;
		return *this;
	}
	Vector4 operator- (const Vector4& v) const
	{
		return Vector4(x-v.x,y-v.y,z-v.z,f-v.f);
	}
	Vector4& operator-= (const Vector4& v)
	{
		x-=v.x;
		y-=v.y;
		z-=v.z;
		f-=v.f;
		return *this;
	}
	Vector4 operator* (const float& n) const
	{
		return Vector4(x*n,y*n,z*n,f*n);
	}
	Vector4 operator/ (const float& n) const
	{
		return Vector4(x/n,y/n,z/n,f/n);
	}
	float operator* (const Vector4& v) const		//Scalar product
	{
		return x*v.x+y*v.y+z*v.z+f*v.f;
	}
};

float determ2(const float& x1, const float& y1, const float& x2, const float& y2)
{
	return x1*y2-y1*x2;
}

float determ3(float x1,float y1,float z1,float x2,float y2,float z2,float x3,float y3,float z3)
{
	return x1*determ2(y2,z2,y3,z3)-y1*determ2(x2,z2,x3,z3)+z1*determ2(x2,y2,x3,y3);
}
float determ3(const Vector3& v1, const Vector3& v2, const Vector3& v3)
{
	return v1.x*determ2(v2.y,v2.z,v3.y,v3.z)-v1.y*determ2(v2.x,v2.z,v3.x,v3.z)+v1.z*determ2(v2.x,v2.y,v3.x,v3.y);
}

//Vector product
Vector3 vec_pr(const Vector3& v1, const Vector3& v2)
{
	Vector3 res;
	res.x=determ2(v1.y,v1.z,v2.y,v2.z);
	res.y=determ2(v1.z,v1.x,v2.z,v2.x);
	res.z=determ2(v1.x,v1.y,v2.x,v2.y);
	return res;
}

//Normal to the plane
Vector3 normal(const point3& p, const point3& p1, const point3& p2)
{
	return vec_pr(p1-p,p2-p);
}

//Vector's module
float module(const Vector3& v)
{
	return float(sqrt(v*v));
}

//Normal to the hyperplane in 4D
Vector4 normal4d(const point4& p, const point4& p1, const point4& p2, const point4& p3)
{
	Vector4 res;
	res.x=determ3(p1.y-p.y,p1.z-p.z,p1.f-p.f,p2.y-p.y,p2.z-p.z,p2.f-p.f,p3.y-p.y,p3.z-p.z,p3.f-p.f);
	res.y=-determ3(p1.x-p.x,p1.z-p.z,p1.f-p.f,p2.x-p.x,p2.z-p.z,p2.f-p.f,p3.x-p.x,p3.z-p.z,p3.f-p.f);
	res.z=determ3(p1.x-p.x,p1.y-p.y,p1.f-p.f,p2.x-p.x,p2.y-p.y,p2.f-p.f,p3.x-p.x,p3.y-p.y,p3.f-p.f);
	res.f=-determ3(p1.x-p.x,p1.y-p.y,p1.z-p.z,p2.x-p.x,p2.y-p.y,p2.z-p.z,p3.x-p.x,p3.y-p.y,p3.z-p.z);
	return res;
}
Vector4 normal4d(const Vector4& a, const Vector4& b, const Vector4& c)
{
	Vector4 res;
	res.x=determ3(a.y,a.z,a.f,b.y,b.z,b.f,c.y,c.z,c.f);
	res.y=-determ3(a.x,a.z,a.f,b.x,b.z,b.f,c.x,c.z,c.f);
	res.z=determ3(a.x,a.y,a.f,b.x,b.y,b.f,c.x,c.y,c.f);
	res.f=-determ3(a.x,a.y,a.z,b.x,b.y,b.z,c.x,c.y,c.z);
	return res;
}

//Vector's module in 4D
float module4d(const Vector4& v)
{
	return float(sqrt(v*v));
}

float cos_angle2(float x1,float y1,float x2,float y2,float x3,float y3)
{
	return (dist2sq(x1,y1,x2,y2)+dist2sq(x3,y3,x2,y2)-dist2sq(x1,y1,x3,y3))/
		2/float(sqrt(dist2sq(x1,y1,x2,y2)*dist2sq(x3,y3,x2,y2)));
}

float cos_angle3(float x1,float y1,float f1,float x2,float y2,float f2,float x3,float y3,float f3)
{
	return (dist3sq(x1,y1,f1,x2,y2,f2)+dist3sq(x3,y3,f3,x2,y2,f2)-dist3sq(x1,y1,f1,x3,y3,f3))/
		2/float(sqrt(dist3sq(x1,y1,f1,x2,y2,f2)*dist3sq(x3,y3,f3,x2,y2,f2)));
}
float cos_angle3(const point3& p1, const point3& p2, const point3& p3)
{
	return (dist3sq(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z)+dist3sq(p3.x,p3.y,p3.z,p2.x,p2.y,p2.z)-
		dist3sq(p1.x,p1.y,p1.z,p3.x,p3.y,p3.z))/
		2/float(sqrt(dist3sq(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z)*dist3sq(p3.x,p3.y,p3.z,p2.x,p2.y,p2.z)));
}

float triangle_perimeter(float x1,float y1,float z1,float x2,float y2,float z2,float x3,float y3,float z3)
{
	double a=sqrt(dist3sq(x1,y1,z1,x2,y2,z2));
	double b=sqrt(dist3sq(x3,y3,z3,x2,y2,z2));
	double c=sqrt(dist3sq(x1,y1,z1,x3,y3,z3));
	return float(a+b+c);
}
float triangle_perimeter(const point3& p1,const point3& p2,const point3& p3)
{
	double a=sqrt(dist3sq(p1,p2));
	double b=sqrt(dist3sq(p3,p2));
	double c=sqrt(dist3sq(p1,p3));
	return float(a+b+c);
}

float triangle_square(float x1,float y1,float z1,float x2,float y2,float z2,float x3,float y3,float z3)
{
	double a=sqrt (dist3sq(x1,y1,z1,x2,y2,z2));
	double b=sqrt (dist3sq(x3,y3,z3,x2,y2,z2));
	double c=sqrt (dist3sq(x1,y1,z1,x3,y3,z3));
	double p=(a+b+c)/2;
	return float(sqrt (p* (p-a) * (p-b) * (p-c) ));
}
float triangle_square(const point3& p1,const point3& p2,const point3& p3)
{
	double a=sqrt(dist3sq(p1,p2));
	double b=sqrt(dist3sq(p3,p2));
	double c=sqrt(dist3sq(p1,p3));
	double p=(a+b+c)/2;
	return float(sqrt (p* (p-a) * (p-b) * (p-c) ));
}
/*float triangle_square(float x1,float y1,float z1,float x2,float y2,float z2,float x3,float y3,float z3)
{
	return module(vec_pr(Vector3(x2-x1,y2-y1,z2-z1),Vector3(x3-x1,y3-y1,z3-z1)))/2;
}
float triangle_square(const point3& p1,const point3& p2,const point3& p3)
{
	return module(vec_pr(Vector3(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z),Vector3(p3.x-p1.x,p3.y-p1.y,p3.z-p1.z)))/2;
}*/

float sin_angle2(float x1,float y1,float x2,float y2,float x3,float y3)
{
	return 2*triangle_square(x1,y1,0.f,x2,y2,0.f,x3,y3,0.f)/
		float(sqrt (dist2sq(x1,y1,x2,y2)*dist2sq(x3,y3,x2,y2)));
}

//Distance between a point and a line in 3D
float dist_point_line(float x,float y,float z,float x1,float y1,float z1,float x2,float y2,float z2)
{
	return 2*triangle_square(x,y,z,x1,y1,z1,x2,y2,z2)/float(sqrt (dist3sq(x1,y1,z1,x2,y2,z2)));
}
float dist_point_line(const point3& p,const point3& p1,const point3& p2)
{
	return 2*triangle_square(p,p1,p2)/float(sqrt(dist3sq(p1,p2)));
}

//Distance between a point and a plane p1p2p3 in 4D
float dist_point_plane4d(const point4& p,const point4& p1,const point4& p2,const point4& p3)
{
	Vector4 norm,norm1;
	Vector4 a,b,c;

	a=Vector4(p1,p2);
	b=Vector4(p1,p3);
	c=Vector4(p1,p);

	norm1=normal4d(a,b,c);
	norm=normal4d(a,b,norm1);

	return float(fabs(norm*c/module4d(norm)));
}

//Distance between a point and a plane
float dist_point_plane(const point3& p,const point3& p1,const point3& p2,const point3& p3)
{
	Vector3 norm=normal(p1,p2,p3);
	return float(fabs(norm*Vector3(p1,p)))/module(norm);
}
float dist_point_plane2(const point3& p,const point3& p1,const point3& p2,const point3& p3)
{
	return float(fabs(determ3(Vector3(p1,p),Vector3(p1,p2),Vector3(p1,p3))))/triangle_square(p1,p2,p3)/2;
}

//Distance between a point and a hyperplane in 4D
float dist_point_hyperplane4d(const point4& p,const point4& p1,const point4& p2,const point4& p3,const point4& p4)
{
	Vector4 norm;
	norm=normal4d(p1,p2,p3,p4);

	return norm*Vector4(p1,p)/module4d(norm);
}

//Verification of (x,y) to be between rays 21 and 23 (from 21 to 23 anticlockwise)
//Dobav' Eps!
int p_in_angle(float x,float y,float x1,float y1,float x2,float y2,float x3,float y3)
{
	float a,b,c;
	double l21,l20,l23;
	int aa=0,bb=0,cc=0;
	const float d=0.01f;

	a=determ2(x1-x2,y1-y2,x-x2,y-y2);
	b=determ2(x-x2,y-y2,x3-x2,y3-y2);
	c=determ2(x1-x2,y1-y2,x3-x2,y3-y2);
	l21=sqrt (dist2sq(x1,y1,x2,y2));
	l20=sqrt (dist2sq(x,y,x2,y2));
	l23=sqrt (dist2sq(x3,y3,x2,y2));

	if ( (fabs(a)/l21/l20<d) && (cos_angle2(x1,y1,x2,y2,x,y)>0.) ) return 0;
	if ( (fabs(b)/l20/l23<d) && (cos_angle2(x3,y3,x2,y2,x,y)>0.) ) return 0;
	if ( (c==0.) && (cos_angle2(x3,y3,x2,y2,x1,y1)>0.) )
	{
		printf ("Error in p_in_angle()\n");
		return -1;
	}

	aa=(a>=0.) ? 1 : 0;
	bb=(b>=0.) ? 1 : 0;
	cc=(c>=0.) ? 1 : 0;
	return ( (aa && bb) || (!cc && ( (!aa && bb) || (aa && !bb) ) ) );
}

//Verification of the projection of a point at the triangle is inside of it in 3D
int pr_in_tr(const point3& p, const point3& p1, const point3& p2, const point3& p3)
{
	Vector3 norm=vec_pr(Vector3(p1,p2),Vector3(p1,p3));
	float E=Eps((fabs(p.x)+fabs(p.y)+fabs(p.z)+fabs(p1.x)+fabs(p1.y)+fabs(p1.z)+
		fabs(p2.x)+fabs(p2.y)+fabs(p2.z)+fabs(p3.x)+fabs(p3.y)+fabs(p3.z))/12);

	Vector3 vec=vec_pr(norm,Vector3(p1,p2));
	if (vec*Vector3(p1,p)<-E) return 0;

	vec=vec_pr(norm,Vector3(p2,p3));
	if (vec*Vector3(p2,p)<-E) return 0;

	vec=vec_pr(norm,Vector3(p3,p1));
	if (vec*Vector3(p3,p)<-E) return 0;

	return 1;
}

//Projection of a point at the plane
point3 pr_p_pl(const point3& p, const point3& p1, const point3& p2, const point3& p3)
{
	Vector3 norm=normal(p1,p2,p3);
	float t=(norm*Vector3(p,p1))/(norm*norm);
	return p+norm*t;
}
point3 pr_p_pl(const point3& p, float A, float B, float C, float D)
{
	Vector3 norm(A,B,C);
	float t=(-D-(norm*Vector3(p)))/(norm*norm);
	return p+norm*t;
}

//Proverka togo, chto proekcia tochki p na pryamuyu zadavaemuyu tochkami p1 i p2, lezhit na otrezke p1p2 (ili v maloi okrestnosti ego kontsov)
//Esli belong=1, znachit izvestno, chto p lezhit na etoi pryamoi
int p_in_segm(const point3& p, const point3& p1, const point3& p2, int belong=0)
{
	float E=Eps((fabs(p.x)+fabs(p.y)+fabs(p.z)+fabs(p1.x)+fabs(p1.y)+fabs(p1.z)+
		fabs(p2.x)+fabs(p2.y)+fabs(p2.z))/9);
	Vector3 v1(p1,p),v2(p2,p);

	if (belong)
	{
		if (v1*v2<=E) return 1;
		else return 0;
	}

	Vector3 vec(p1,p2);
	if (vec*v1>=-E && vec*v2<=E) return 1;
	else return 0;
}

//Proverka togo, chto tochka p lezhit vnutri tetraedra p1p2p3p4 (ili na ego granitse)
//Dobav' Eps!
int p_in_tet(const point3& p, const point3& p1, const point3& p2, const point3& p3, const point3& p4)
{
	point3 v[4];
	Vector3 norm,pp=Vector3(p);

	v[0]=p1;
	v[1]=p2;
	v[2]=p3;
	v[3]=p4;

	for (int i=0; i<4; i++)
	{
		norm=normal(v[i],v[(i+1)%4],v[(i+2)%4]);
		if ((norm*Vector3(v[i],v[(i+3)%4]))*(norm*Vector3(v[i],p))<0) return 0;
	}
	
	return 1;
}

//point of intersection of 2 lines
point2 p_of_inter(float x1,float y1,float x2,float y2,float x3,float y3,float x4,float y4)
{
	point2 pp;
	float d,t;

	t=determ2(x3-x1,y3-y1,x4-x3,y4-y3);
	d=determ2(x2-x1,y2-y1,x4-x3,y4-y3);

	float E=Eps((fabs(x1)+fabs(y1)+fabs(x2)+fabs(y2)+fabs(x3)+fabs(y3)+fabs(x4)+fabs(y4))/8);
	E*=E;

	if (fabs(d)>E)
	{
		pp.x=x1+(x2-x1)*t/d;
		pp.y=y1+(y2-y1)*t/d;
	}
	else
	{
		pp.x=3.3e38f;
		pp.y=3.3e38f;
	}

	return pp;
}

//Point of intersection of 2 lines in 3D
//Pryamye ne dolzhny byt' parallel'ny!
//Pryamye dolzhny lezhat' v odnoi ploskosti, inache poluchaemaya tochka budet lezhat' tol'ko na 1-i pryamoi, 
//no budet smeshchena vdol' odnoi iz koordinatnykh osei otnositel'no nekotoroi tochki 2-i pryamoi
point3 p_of_inter3D(const point3& p1,const point3& p2,const point3& q1,const point3& q2)
{
	//Uravnenie 1-i pryamoi: p=p1+(p2-p1)*t=p1+a*t
	//Uravnenie 2-i pryamoi: q=q1+(q2-q1)*u=q1+b*u
	Vector3 a(p1,p2),b(q1,q2),w(p1,q1);
	float t;

	float d_xy=determ2(a.x,b.x,a.y,b.y);
	float d_xz=determ2(a.x,b.x,a.z,b.z);
	float d_yz=determ2(a.y,b.y,a.z,b.z);

	if (fabs(d_xy)>=fabs(d_xz) && fabs(d_xy)>=fabs(d_yz))
	{
		t=determ2(w.x,b.x,w.y,b.y)/d_xy;
	}
	else if (fabs(d_xz)>=fabs(d_xy) && fabs(d_xz)>=fabs(d_yz))
	{
		t=determ2(w.x,b.x,w.z,b.z)/d_xz;
	}
	else
	{
		t=determ2(w.y,b.y,w.z,b.z)/d_yz;
	}

	return p1+a*t;
}

//Tochka peresecheniya pryamoy i ploskosti
//Pryamaya zadaetsya 2-ya tochkami, ploskost' - koeffitsientami iz uravneniya Ax+By+Cz+D=0
point3 line_plane(float x1,float y1,float z1,float x2,float y2,float z2,float A,float B,float C,float D)
{
	float a=x2-x1;
	float b=y2-y1;
	float c=z2-z1;
	float t=A*a+B*b+C*c;
	point3 result;

	float E=Eps((fabs(x1)+fabs(y1)+fabs(z1)+fabs(x2)+fabs(y2)+fabs(z2))/6);

	if (fabs(t)<=E)
	{
		result.x=3.3e38f;
		result.y=3.3e38f;
		result.z=3.3e38f;

		return result;
	}

	t=-(A*x1+B*y1+C*z1+D)/t;
	result.x=x1+a*t;
	result.y=y1+b*t;
	result.z=z1+c*t;

	return result;
}
//Pryamaya i ploskost' zadayutsya tochkami
point3 line_plane(const point3& p1,const point3& p2,const point3& q1,const point3& q2,const point3& q3)
{
	Vector3 vec(p1,p2);
	Vector3 norm=normal(q1,q2,q3);
	float t=norm*vec;
	point3 result;

	float E=Eps((fabs(p1.x)+fabs(p1.y)+fabs(p1.z)+fabs(p2.x)+fabs(p2.y)+fabs(p2.z)+
		fabs(q1.x)+fabs(q1.y)+fabs(q1.z)+fabs(q2.x)+fabs(q2.y)+fabs(q2.z)+
		fabs(q3.x)+fabs(q3.y)+fabs(q3.z))/15);

	if (fabs(t)<=E)
	{
		result.x=3.3e38f;
		result.y=3.3e38f;
		result.z=3.3e38f;

		return result;
	}

	t=(norm*Vector3(p1,q1))/t;
	return p1+vec*t;
}

//Ob'em tetraedra
float tetrahedron_volume(float x1,float y1,float z1,float x2,float y2,float z2,float x3,float y3,float z3,float x4,float y4,float z4)
{
	return float(fabs(determ3(x2-x1,y2-y1,z2-z1,x3-x1,y3-y1,z3-z1,x4-x1,y4-y1,z4-z1)));
}

//Maksimumy i minimumy
float maxim(float x,float y)
{
	return x>y ? x : y;
}
float maxim(float x,float y,float z)
{
	return maxim(x,maxim(y,z));
}
int maxim(int a,int b)
{
	return a>b ? a : b;
}

float minim(float x,float y)
{
	return x<y ? x : y;
}
int minim(int a,int b)
{
	return a<b ? a : b;
}

//Rasstoyanie mezhdu tochkoi i otrezkom
float dist_point_segm(const point3& p,const point3& p1,const point3& p2)
{
	float h=3.3e38f,temp;
	float E=Eps((fabs(p.x)+fabs(p.y)+fabs(p.z)+fabs(p1.x)+fabs(p1.y)+fabs(p1.z)+
		fabs(p2.x)+fabs(p2.y)+fabs(p2.z))/9);

	if (p_in_segm(p,p1,p2))
	{
		h=dist_point_line(p,p1,p2);
		if (h<=E) return 0;
		return h;
	}

	temp=float(sqrt(dist3sq(p,p1)));
	if (temp<=E) return 0;
	h=float(sqrt(dist3sq(p,p2)));
	if (h<=E) return 0;

	if (h>temp) h=temp;

	return h;
}

//Rasstoyanie mezhdu 2-ya otrezkami
float dist_segm_segm(const point3& p1,const point3& p2,const point3& q1,const point3& q2)
{
	float H=3.3e38f,temp;
	Vector3 v1(p1,p2),v2(q1,q2);
	float E=Eps((fabs(p1.x)+fabs(p1.y)+fabs(p1.z)+fabs(p2.x)+fabs(p2.y)+fabs(p2.z)+
		fabs(q1.x)+fabs(q1.y)+fabs(q1.z)+fabs(q2.x)+fabs(q2.y)+fabs(q2.z))/12);
	int i,j;
	point3 points[4];
	points[0]=p1;
	points[1]=p2;
	points[2]=q1;
	points[3]=q2;

	//Proveryaem parallel'nost' pryamykh
	Vector3 norm=vec_pr(v1,v2);
	if (module(norm)<=E)
	{
		if (p_in_segm(q1,p1,p2) || p_in_segm(q2,p1,p2))
		{
			H=dist_point_line(q1,p1,p2);
			if (H<=E) return 0;
			return H;
		}
		else
		{
			for (i=0; i<2; i++)
			{
				for (j=2; j<4; j++)
				{
					temp=float(sqrt(dist3sq(points[i],points[j])));
					if (H>temp) H=temp;
					if (H<=E) return 0;
				}
			}
			return H;
		}
	}

	//Pryamye neparallel'ny
	float D1=-(norm*Vector3(p1));
	float D2=-(norm*Vector3(q1));
	point3 pr1=pr_p_pl(q1,norm.x,norm.y,norm.z,D1);
	point3 pr2=pr_p_pl(q2,norm.x,norm.y,norm.z,D1);
	point3 pp=p_of_inter3D(pr1,pr2,p1,p2);
	if (p_in_segm(pp,pr1,pr2,1) && p_in_segm(pp,p1,p2,1))
	{
		H=float(fabs(D2-D1)/module(norm));
		if (H<=E) return 0;
		return H;
	}

	temp=dist_point_segm(p1,q1,q2);
	if (H>temp) H=temp;
	if (H<=E) return 0;
	temp=dist_point_segm(p2,q1,q2);
	if (H>temp) H=temp;
	if (H<=E) return 0;
	temp=dist_point_segm(q1,p1,p2);
	if (H>temp) H=temp;
	if (H<=E) return 0;
	temp=dist_point_segm(q2,p1,p2);
	if (H>temp) H=temp;
	if (H<=E) return 0;

	return H;
}

//Rasstoyanie mezhdu otrezkom i treugol'nikom
float dist_segm_tr(const point3& p1,const point3& p2,const point3& q1,const point3& q2,const point3& q3)
{
	int i,j;
	float A,B,C,D,D1,d,temp1,temp2;
	point3 q[3],pp;
	q[0]=q1;
	q[1]=q2;
	q[2]=q3;

	Vector3 norm=normal(q1,q2,q3);

	A=norm.x;
	B=norm.y;
	C=norm.z;
	D=-(norm*Vector3(q1));

	double mmm=fabs(p1.x)+fabs(p1.y)+fabs(p1.z)+fabs(p2.x)+fabs(p2.y)+fabs(p2.z);
	for (i=0; i<3; i++)
	{
		mmm+=fabs(q[i].x)+fabs(q[i].y)+fabs(q[i].z);
	}
	float E=Eps(mmm/15);

	Vector3 dir_vec(p1,p2);		//Napravlyayuschii vektor pryamoi otrezka

	//Proveryaem parallel'nost' pryamoi otrezka i ploskosti
	int is_par=0;
	if (fabs(norm*dir_vec)<=E)
	{
		is_par=1;

		//pr1-pr2 - proektsiya otrezka na ploskost'
		point3 pr1=pr_p_pl(p1,q1,q2,q3);
		point3 pr2=pr_p_pl(p2,q1,q2,q3);

		//Proveryaem peresechenie proektsii s rebrami treugol'nika
		int inter=0;	//Peresekaet li proektsiya treugol'nik
		for (i=0; i<2 && inter==0; i++)
		{
			for (j=i+1; j<3; j++)
			{
				//Proveryaem parallel'nost' pryamoi pr1-pr2 i rebra treugol'nika
				if (module(vec_pr(Vector3(pr1,pr2),Vector3(q[i],q[j])))<=E)
				{
					//Proveryaem sovpadenie tochek
					if (pr1==q[i] || pr1==q[j] || pr2==q[i] || pr2==q[j])
					{
						inter=1;
						break;
					}

					//Proveryaem sovpadenie pryamykh i nalozhenie otrezkov
					if (module(vec_pr(Vector3(q[i],pr1),Vector3(q[i],q[j])))<=E && 
						(p_in_segm(pr1,q[i],q[j],1) || p_in_segm(pr2,q[i],q[j],1)))
					{
						inter=1;
						break;
					}
				}
				else
				{
					pp=p_of_inter3D(pr1,pr2,q[i],q[j]);		//Tochka peresecheniya pryamykh
					if (p_in_segm(pp,pr1,pr2,1) && p_in_segm(pp,q[i],q[j],1))
					{
						inter=1;
						break;
					}
				}
			}
		}

		if (inter)
		{
			D1=-(norm*Vector3(p1));
			d=float(fabs(D1-D)/sqrt(A*A+B*B+C*C));
			if (d<=E) return 0;
			return d;
		}
	}

	d=3.3e38f;

	if (!is_par)
	{
		//Nakhodim tochku peresecheniya pryamoi otrezka i ploskosti treugol'nika
		pp=line_plane(p1,p2,q1,q2,q3);
		//Proveryaem prinadlezhnost' pp otrezku i treugol'niku
		if (p_in_segm(pp,p1,p2,1) && pr_in_tr(pp,q1,q2,q3))
		{
			return 0;
		}
		else
		{
			//Nakhodim konets otrezka, kotoryi yavlyaetsya blizhaishim k ploskosti treugol'nika, i 
			//proektsiya kotorogo popadaet v treugol'nik
			temp1=dist_point_plane(p1,q1,q2,q3);
			temp2=dist_point_plane(p2,q1,q2,q3);
			if (temp1<=temp2)
			{
				if (pr_in_tr(p1,q1,q2,q3)) d=temp1;
			}
			else
			{
				if (pr_in_tr(p2,q1,q2,q3)) d=temp2;
			}
		}
	}

	if (d<=E) return 0;

	//Nakhodim minimum rasstoyanii mezhdu otrezkom i rebrami treugol'nika
	for (i=0; i<2; i++)
	{
		for (j=i+1; j<3; j++)
		{
			temp1=dist_segm_segm(p1,p2,q[i],q[j]);
			if (d>temp1) d=temp1;
			if (d==0) return 0;
		}
	}

	if (d<=E) return 0;

	return d;
}

//Rasstoyanie mezhdu 2-ya treugol'nikami
float dist_tr_tr(float** tr1, float** tr2)
{
	int i,j;
	float D,temp;
	point3 p[3],q[3];

	for (i=0; i<3; i++)
	{
		p[i].x=tr1[i][0];
		p[i].y=tr1[i][1];
		p[i].z=tr1[i][2];

		q[i].x=tr2[i][0];
		q[i].y=tr2[i][1];
		q[i].z=tr2[i][2];
	}

	Vector3 norm1=normal(p[0],p[1],p[2]);
	Vector3 norm2=normal(q[0],q[1],q[2]);

	float mm=0;
	for (i=0; i<3; i++)
	{
		for (j=0; j<3; j++)
		{
			mm+=float(fabs(tr1[i][j]));
		}
	}
	for (i=0; i<3; i++)
	{
		for (j=0; j<3; j++)
		{
			mm+=float(fabs(tr2[i][j]));
		}
	}

	float E=Eps(mm/18);

	D=3.3e38f;
	for (i=0; i<2; i++)
	{
		for (j=i+1; j<3; j++)
		{
			if (D>(temp=dist_segm_tr(p[i],p[j],q[0],q[1],q[2])))
			{
				D=temp;
				if (D==0) return 0;
			}
			if (D>(temp=dist_segm_tr(q[i],q[j],p[0],p[1],p[2])))
			{
				D=temp;
				if (D==0) return 0;
			}
		}
	}

	if (D<=E) return 0;

	return D;
}

template<class T> class Matrix2D
{
	T** matr;
	int n,m;
public:
	Matrix2D(int nn, int mm)
	{
		matr=new T*[nn];
		for (int i=0; i<nn; i++)
		{
			matr[i]=new T[mm];
		}
		n=nn;
		m=mm;
	}

	Matrix2D(const Matrix2D& mtr)
	{
		n=mtr.n;
		m=mtr.m;
		matr=new T*[n];
		for (int i=0; i<n; i++)
		{
			matr[i]=new T[m];
			for (int j=0; j<m; j++)
			{
				matr[i][j]=mtr.matr[i][j];
			}
		}
	}

	~Matrix2D()
	{
		for (int i=0; i<n; i++)
		{
			delete[] matr[i];
		}
		delete[] matr;
	}

	Matrix2D& operator= (const Matrix2D& mtr)
	{
		if (n!=mtr.n || m!=mtr.m)
		{
			for (int i=0; i<n; i++)
			{
				delete[] matr[i];
			}
			delete[] matr;
			matr[0][0]=mtr.matr[0][0];		//Exit
		}
		else
		{
			for (int i=0; i<n; i++)
			{
				for (int j=0; j<m; j++)
				{
					matr[i][j]=mtr.matr[i][j];
				}
			}
		}

		return *this;
	}

	Matrix2D operator+ (const Matrix2D& mtr) const
	{
		Matrix2D res(n,m);
		if (n!=mtr.n || m!=mtr.m)
		{
			for (int i=0; i<n; i++)
			{
				delete[] matr[i];
			}
			delete[] matr;
			matr[0][0]=mtr.matr[0][0];		//Exit with error
		}
		else
		{
			for (int i=0; i<n; i++)
			{
				for (int j=0; j<m; j++)
				{
					res.matr[i][j]=matr[i][j]+mtr.matr[i][j];
				}
			}
		}

		return res;
	}

	Matrix2D& operator+= (const Matrix2D& mtr)
	{
		if (n!=mtr.n || m!=mtr.m)
		{
			for (int i=0; i<n; i++)
			{
				delete[] matr[i];
			}
			delete[] matr;
			matr[0][0]=mtr.matr[0][0];		//Exit
		}
		else
		{
			for (int i=0; i<n; i++)
			{
				for (int j=0; j<m; j++)
				{
					matr[i][j]+=mtr.matr[i][j];
				}
			}
		}

		return *this;
	}

	Matrix2D operator* (const Matrix2D& mtr) const
	{
		Matrix2D res(n,mtr.m);
		if (m!=mtr.n)
		{
			for (int i=0; i<n; i++)
			{
				delete[] matr[i];
			}
			delete[] matr;
			matr[0][0]=mtr.matr[0][0];		//Exit with error
		}
		else
		{
			for (int i=0; i<n; i++)
			{
				for (int j=0; j<mtr.m; j++)
				{
					res.matr[i][j]=0;
					for (int ii=0; ii<m; ii++)
					{
						res.matr[i][j]+=matr[i][ii]*mtr.matr[ii][j];
					}
				}
			}
		}

		return res;
	}

	T*& operator[] (const int& n)
	{
		return matr[n];
	}

	int lin () const
	{
		return n;
	}

	int col () const
	{
		return m;
	}

	Matrix2D& exch (int i, int j)		//Pomenyat' v matritse stroki i i j
	{
		T temp;
		if (i!=j)
		{
			for (int ii=0; ii<m; ii++)
			{
				temp=matr[i][ii];
				matr[i][ii]=matr[j][ii];
				matr[j][ii]=temp;
			}
		}

		return *this;
	}

	Matrix2D inv () const
	{
		int i,j,k;
		Matrix2D err(0,0);

		if (n!=m)
		{
			return err;
		}

		//Prislonyaem sprava edinichnuyu matritsu
		Matrix2D mtr(n,2*n);
		for (i=0; i<n; i++)
		{
			for (j=0; j<n; j++)
			{
				mtr.matr[i][j]=matr[i][j];
			}
			for (j=0; j<n; j++)
			{
				mtr.matr[i][j+n]=(i==j?T(1.):0);
			}
		}

		double mm=0;
		for (i=0; i<n; i++)
		{
			for (j=0; j<n; j++)
			{
				mm+=fabs(matr[i][j]);
			}
		}
		float E=Eps(mm/n/n);

		int cur_line=0;
		T temp,max;
		while (cur_line<n)
		{
			//Ishem stroku k s maksimalnym po modulyu elementom mtr.matr[k][cur_line]
			k=cur_line;
			max=T(fabs(mtr.matr[k][cur_line]));
			i=cur_line+1;
			while (i<n)
			{
				temp=T(fabs(mtr.matr[i][cur_line]));
				if (temp>max)
				{
					k=i;
					max=temp;
				}
				i++;
			}

			if (max<=E)
			{
				return err;
			}

			//Menyaem mestami stroki k i cur_line
			if (k!=cur_line)
			{
				for (int ii=0; ii<2*n; ii++)
				{
					temp=mtr.matr[k][ii];
					mtr.matr[k][ii]=mtr.matr[cur_line][ii];
					mtr.matr[cur_line][ii]=temp;
				}
			}

			//Izmenyaem stroku tak, chtoby element mtr.matr[cur_line][cur_line] stal edinitsei
			temp=mtr.matr[cur_line][cur_line];
			mtr.matr[cur_line][cur_line]=T(1.);
			for (i=cur_line+1; i<2*n; i++)
			{
				mtr.matr[cur_line][i]/=temp;
			}

			//V ostalnykh strokakh stolbtsa cur_line nuzhno poluchit' nuli
			for (i=0; i<n; i++)
			{
				if (i==cur_line)
				{
					continue;
				}
				temp=mtr.matr[i][cur_line];
				mtr.matr[i][cur_line]=0;
				for (j=cur_line+1; j<2*n; j++)
				{
					mtr.matr[i][j]-=mtr.matr[cur_line][j]*temp;
				}
			}

			cur_line++;
		}

		Matrix2D res(n,n);
		for (i=0; i<n; i++)
		{
			for (j=0; j<n; j++)
			{
				res.matr[i][j]=mtr.matr[i][j+n];
			}
		}

		return res;
	}
};

#endif

/*
//Poisk determinanta i obratnoi matritsy W dlya matritsy S 4x4 (metod imeet ochen' plohuyu shodimost', kogda dopolnitel'nye minory <<1)
determ=0;
for (i=0; i<4; i++)
{
	determ+=S[0][i]*determ3(S[1][(i+1)%4],S[1][(i+2)%4],S[1][(i+3)%4],S[2][(i+1)%4],
		S[2][(i+2)%4],S[2][(i+3)%4],S[3][(i+1)%4],S[3][(i+2)%4],S[3][(i+3)%4]);
}
for (i=0; i<4; i++)
{
	for (j=0; j<4; j++)
	{
		W[i][j]=((i+j)%2?-1:1)*determ3(S[(j+1)%4][(i+1)%4],S[(j+1)%4][(i+2)%4],S[(j+1)%4][(i+3)%4],
				S[(j+2)%4][(i+1)%4],S[(j+2)%4][(i+2)%4],S[(j+2)%4][(i+3)%4],
				S[(j+3)%4][(i+1)%4],S[(j+3)%4][(i+2)%4],S[(j+3)%4][(i+3)%4])/determ;
	}
}
*/
